Thermalization of a pump-excited Mott insulator 
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We use nonequilibrium dynamical mean-field theory in combination with a recently implemented 
strong-coupling impurity solver to investigate the relaxation of a Mott insulator after a laser excita- 
tion with frequency comparable to the Hubbard gap. The time evolution of the double occupancy 
exhibits a crossover from a strongly damped transient at short times towards an exponential ther- 
malization at long times. In the limit of strong interactions, the thermalization time is consistent 
with the exponentially small decay rate for artificially created doublons, which was measured in ul- 
tracold atomic gases. When the interaction is comparable to the bandwidth, on the other hand, the 
double occupancy thermalizes within a few times the inverse bandwidth along a rapid thermaliza- 
tion path in which the exponential tail is absent. Similar behavior can be observed in time-resolved 
photoemission spectroscopy. Our results show that a simple quasi-equilibrium description of the 
electronic state breaks down for pump-excited Mott insulators characterized by strong interactions. 



(N 



i 

h 



i 

C 

o 
o 



> 
in 
q 

o 



X 



I. INTRODUCTION 

State-of-the-art pump-probe experiments use ultra- 
short light pulses to control and probe the complex 
interplay between electronic, magnetic and lattice de- 
grees of freedom in strongly correlated materials The 
Mott metal-insulator transition, which is the paradigm 
for phase transitions in correlated electron systems, has 
been investigated in several time-resolved experiments. 
For example, a photo-induced metal-insulator transition 
and the ultrafast recovery of the Mott gap was ob- 
served in a charge-transfer insulator by using reflectiv- 
ity measurements p and in lT-TaS2 using time-resolved 
photoemission spectroscopy^ Recently, Wall et al. suc- 
ceeded in measuring electronic dynamics in an organic 
charge-transfer insulator on the timescalc of the inverse 
bandwidth h/W A 

Photo-excited phases in materials with several com- 
peting interactions have been reached along nonthermal 
pathways^ and their properties can be quite different 
from those of any known equilibrium state of the same 
system^^ On the other hand, in many experiments the 
electronic state has been interpreted in terms of the two- 
temperature model^ whose basic assumption is that elec- 
trons thermalize rapidly (on the timescale of the inverse 
bandwidth) and thereafter can be assigned a tempera- 
ture which is different from the lattice temperature. De- 
viations from the two-temperature behavior have been 
observed long ago in noble metals, where the scattering 
between quasiparticles is weak and thermalization times 
can range up to a few hundred femtoseconds £ In strongly 
interacting quantum systems, thermalization is theoreti- 
cally not well understood, and it is thus unclear whether 
a homogeneous electron system can be stuck in a non- 
thermal state on timescales much longer than h/W . 

Nonequilibrium experiments on ultracold atomic 
gases^^ have recently stimulated considerable theo- 
retical effort to resolve some of the questions related 
to the thermalization of isolated quantum-many body 
systemsAi In several models^— it was found that the 



relaxation after a sudden quench of the Hamiltonian sen- 
sitively depends on the interaction parameter. In the 
fermionic Hubbard model, e.g., rapid thermalization of 
the momentum distribution (within t ~ h/W) occurs 
after quenches from the noninteracting ground state to 
a very narrow interaction regime at intermediate cou- 
pling (U ~ 0.8M 7 for the semi-elliptic density of states)^ 
while after quenches to weak and strong coupling one 
observes the formation of a quasistationary nonthermal 
state which persists for t 3> H/W i 16 The latter phe- 
nomenon, which is called prethermalization, can be re- 
lated to the non-ergodic behavior of integrable models^ 
because in both cases the dynamics is constrained by ei- 
ther exact or approximate constants of motion*^ 

In this paper we turn from interaction quenches to the 
pump-probe setup and study the relaxation of a Mott 
insulator after a short laser excitation. The analysis is 
performed for the Hubbard model, 



y,cr=t,4- i 

which describes fermionic particles that can hop between 
the sites of a crystal lattice (with hopping amplitude ty) 
and interact with each other through a local Coulomb re- 
pulsion U. In the phase diagram of the Hubbard model 
there is a metal-insulator transition at low temperatures 
and half- filling (n-j- = = ^). For the present investiga- 
tion we focus on higher temperatures, where this phase 
transition turns into a smooth crossover, and study the 
relaxation after a short pump-pulse that excites carriers 
over the Mott gap. Although the equilibrium states be- 
fore the excitation are then continuously connected, the 
response to the pulse turns out to be strongly dependent 
on the value of U, 

For all values of U, the system approaches a thermal 
state, which is at elevated temperature because the en- 
ergy is increased during the pump and conserved there- 
after. However, just as for the interaction quench« the 
thermalization time r t h depends sensitively on U . Al- 
ready for moderately strong interactions, where the Hub- 
bard gap is still comparable to the bandwidth, rth is much 
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larger than h/W, which puts a question mark behind the 
applicability of a simple two-temperature model in this 
case. At intermediate coupling, on the other hand, ther- 
malization indeed happens on the timescale of the inverse 
bandwidth, but it then follows a qualitatively different 
path which may be referred to as "rapid thermalization" . 

To study the dynamics of the Hubbard model, we 
use the dynamical mean-field theory^ which provides 
a mapping of the lattice model onto a single im- 
purity model and is exact in the limit of infinite 
dimensions^ The reformulation of DMFT within the 
Keldysh framewor k 21 ' 22 allows to apply this approach to 
nonequilibrium problems in which a system is initially 
prepared in a thermal equilibrium state and subsequently 
exposed to some perturbation. The biggest challenge for 
the method is the solution of a quantum impurity model 
under quite general nonequilibrium conditions. A gener- 
alization of the auxiliary-field weak-coupling continuous- 
time quantum Monte Carlo method 2 ^ has been used 
successfully to study the short-time dynamics of the 
Hubbard model ) 14 ' 15 ' 24 ' 25 but the dynamical sign prob- 
lem prevents an investigation of the long-time behav- 
ior, in particular for large values of U which arc of in- 
terest for the current investigation of Mott-insulating 
states. Possibilities to avoid the difficulties of real- 
time Monte Carlo methods include the superperturba- 
tion theory^ and perturbative approaches! 15 ' 27 For the 
present study we use the self-consistent strong-coupling 
expansion^ which is a generalization of the non-crossing 
approximation 2 ^ (NCA) to higher orders and to the 
Keldysh contour. It works particularly well in the Mott 
insulating phase, and the comparison of results from dif- 
ferent orders (up to third order) provides at least some 
estimate of the error. 

The outline of this paper is as follows: In Section 
Hi! we give some details of the DMFT framework used 
here, in particular the solution of the lattice Dyson equa- 
tion in real (Keldysh) time. We then present results for 
the relaxation of the double occupancy after a pump- 
excitation of a Mott insulator (Sec. Illlj) and identify the 
corresponding signatures in time-resolved photoemission 
spectroscopy (Sect. H*V|) . A conclusion is given at the end 
(Sec. E| . 



II. THE MODEL 

A. Hubbard Model in an external electric field 

In this paper we investigate the Hubbard model 
with an additional spatially homogeneous, but time- 
dependent electric field, 



E(t) =nE Q sm(nt)<f>(t), 



(2) 



corresponding to a few-cycle pump-pulse with amplitude 
Eq, polarization h, frequency f2, and a pulse envelope 
4>(t). (The assumption of spatial homogeneity is ap- 
propriate for optical frequencies with wavelength much 



larger than the lattice spacing.) Within a gauge without 
scalar potential, E is determined by the vector poten- 
tial through the relation E = — \dtA, and the latter is 
incorporated into the Hamiltonian ([TJ via a Peierls sub- 
stitution in the hopping amplitudes tij. This leads to a 
shift of the band energies e^, i.e., the hopping part of 
Eq. ([TJ becomes time-dependent, 



H(t) 



hk(t)n ka + Uy^(ni t -±)(n u -±), (3) 



with hk(t) = efc_^yi( t ). In the following we choose a 
hypercubic lattice in the limit of infinite dimensions with 
a Gaussian density of states 
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The unit for energy is given by the variance W of the 
density of states, which will be loosely referred to as the 
bandwidth, time is measured in units of h/W , and the 
unit of the electric field is given by W/ea, where — e is 
the electronic charge and a is the lattice spacing. In 
these units, the critical end-point of the first-order Mott 
transition line in the phase diagram of the paramagnetic 
half-filled Hubbard model, which is obtained quite accu- 
rately at second order in the self-consistent hybridization 
expansion^ is located at U ~ 3.1 and T « 0.02. 

The DMFT equations for this setup have been dis- 
cussed in detail in Ref. [22| for a polarization h = 
(1,1,. ..,1)* along the body diagonal of the cubic unit 
cell, and they will not be repeated in detail. Because 
A(t) enters in the effective impurity action only indirectly 
via the hybridization function, we can apply the self- 
consistent hybridization expansion as an impurity solver 
without modifications^ If not stated otherwise we use 
the second order of that approximation as an impurity 
solver (one-crossing approximation, OCA), which is reli- 
able in the insulating phase and in the crossover regime^ 



B. Solution of the DMFT self-consistency 

Within the DMFT framework, a repeated solution of 
the lattice Dyson equation and the Dyson equation of 
the impurity model is required (for details, see, e.g., 
Refs. [22| and [HI). These equations are in essence integro- 
differential equations on the Keldysh contour C which 
runs from to time i ma x (i-e., the largest time of inter- 
est) along the real axis, back to 0, and finally to —i/3 
along the imaginary time axis^ Since there are in prin- 
ciple several ways to rephrase and solve these equations, 
we use the remaining paragraphs of this section to de- 
scribe a scheme which we found to be reliable. The same 
scheme has been used previously to study the dielectric 
breakdown of the Mott insulator.— 

The following approach is mainly designed to satisfy 
two properties: (i) It is easy to implement a higher order 
scheme in the timcstcp At to reduce the discretization 
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error. Wc found that this considerably increases the ac- 
curacy for a given number of timesteps N on C, which in 
turn strongly reduces the required storage [0(N 2 )] and 
computational effort [0(N 3 ) for NCA, 0(N 4 ) for OCA], 
(ii) The approach exactly satisfies the causal structure of 
the equations, i.e., the solution at a given time is com- 
puted using only information from previous times. This 
makes the solution of the nonequilibrium DMFT equa- 
tions equivalent to a time-propagation that starts from 
the initial equilibrium DMFT solution, rather than an 
iteration of the equations on the full contour C. 

To design such a scheme we rewrite the integro- 
differential equations on C explicitly for the imaginary- 
time, real-time, and mixed imaginary-time /real-time 
components of the two-time propagators G(t, t'). This 
step has been described in detail in Ref. [[a. The re- 
sulting equations are integral equations of Volterra type, 
whose solution can be easily designed in a way to satisfy 
the requirements (i) and (ii) above. We found that two 
types of integral equations are particularly well-behaved 
(both are to be solved for X): The integral equation 

(I - A)*X = B, (5) 

and the integro-differential equation (basically a Dyson 
equation) 

[id t - h - A] * X = I, (6) 

where h(t,t') = 8c(t,t')h(t) is a function of the physical 
time t only. Here and in the following, bold quantities 
X(t,t') denote two-time functions with arguments on C, 
[A * B](t,t') = J c dtA(t,t)B(t,t') is the contour con- 
volution, and I(t,t') = 6c(t,t') is the delta function on 
C. When written as explicit integral equations for the 
real and imaginary-time components, Eq. ([3]) leads to 
Volterra equations of the second kind. In contrast, the 
slightly modified equation A * X = B leads to Volterra 
equations of the first kind, which tend to be unstable for 
high-order algorithms. In the remainder of this section 
we show how to rewrite the DMFT selfconsistency equa- 
tions (appropriate for the NCA-type impurity solver) 
only in terms of equations of type 

Nonequilibrium DMFT corresponds to a self-consistent 
solution of the following set of equations: (i) Solution of 
the impurity problem, which in the case of the strong- 
coupling expansion requires an evaluation of the diagram- 
matic expression for the local Greenfunction in terms of 
the hybridization function A^i 

G = G NCA ,...[A], (7) 

(ii) , the Dyson equation of the impurity model 

[id t + n - £ - A] * G = I, (8) 

(iii) , the lattice Dyson equation with the time-dependent 
band-energies hk(t) from Eq. ^ 

[id t + H - E - h h ] * G k = I, (9) 



and (iv), the momentum sum 

Y, G k = G. (10) 

k 

To reformulate this set of equations only in terms of 
equations of type ([5]) we introduce the quantity Z = 
[id t + — £] _1 , conjugate to the self-energy, and the mo- 
mentum sums Gi = J^k^kGk, G\ = J^k^khk, and 
G2 = J2k hkGkhk in addition to Eq. (|10p . By summing 
Eq. ([5]) over fc (using ^ fc hk = 0) and comparing with 
Eq. © one can show that G\ — A*G, G\ — G * A, and 
that Eqs. (|8j)- (|10p arc equivalent to 

[I + G\]*Z = G, (11) 
[I - F k ] *G k = Z, (12) 
[I + G 1 ]*A = G 2 , (13) 

where Fk{t,t') = Z(t,t')hk(t'). In practice, we deter- 
mine G from A via the self-consistent hybridization ex- 
pansion, perform the convolution G * A to obtain an 
estimate for G\, then solve Eq. (JTT]) for Z, Eq. (fTJ]) for 
Gk on a given set of fe-points, perform the momentum 
sums to obtain G, G\, G\, and G2, and finally solve 
Eq. (flB"!) for A. This cycle is performed consecutively for 
each timestep (starting from an extrapolation of A from 
the previous timestep), which usually does not require 
more than two or three iterations to reach a converged 
solution. 



III. THERM ALIZATION OF THE 
PUMP-EXITED MOTT INSULATOR 

In this section we analyze the dynamics of Mott- 
insulating and crossover states in the half-filled Hubbard 
model after pump-pulses with frequencies of the order 
of the gap. The system is initially prepared in an equi- 
librium state at given temperature T = 1//3, and ex- 
cited by a few-cycle electrical field pulse [Eq. ©] that 
has a Gaussian pulse envelope <j)(t) = exp[— (t — 2) 2 ] and 
is restricted to time < t < 4 (Fig. []Ji). Figure [T|b 
shows the typical time dependence of the total energy 
E{t) = (H(t)) during and after the pump-excitation for 
an initially Mott-insulating state (U = 5, /3 = 5). During 
the pump, £(t) increases with time due to the absorption 
of light, while it is constant or t > 4 because the sys- 
tem is not coupled to external thermal reservoirs. Note 
that energy conservation comes out correct only if all 
approximations that are made for the solution of the im- 
purity model are conserving in the sense of Kadanoff and 
Baym, which is true for the self-consistent hybridization 
expansion. More precisely, energy conservation implies 
the relation £ (t) = E(t)(j(t)) , which can serve as a first 
check for the numerics (j = e ^2 krT rikaVk is the current 
operator, with band velocities Vk — ftr 1 dk^k- ^- a) ■ 

It is the basic assumption of the two-tcmpcraturc 
model that due to the electron-electron interaction elec- 
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FIG. 1: Pump-excitation of a Mott insulator at /3 = 5, U = 5. 

a) The electric field used for the data of this plot. The two 
frequencies fl = 5n/4 ~ 3.927 and = 57r/2 « 7.854 are com- 
parable to the Hubbard gap and larger than the Hubbard gap, 
respectively, and the amplitudes Eo = 0.9679 (SI = 5n/4) and 
-Bo = 3.6701 (Q = 5n/2) are chosen such that T e s = 0.5, i.e., 
the energy of the pump-excited system for t > 4 is the same as 
in a system in thermal equilibrium at temperature T = T e g. 

b) Time-dependence of the total energy £ . c) Time depen- 
dence of the double occupancy d(t). The horizontal dotted 
line indicates the double occupancy d(T e ff) in an equilibrium 
state at temperature T = T e g. 



trons in a solid thermalize so fast that the lattice dy- 
namics is not important for that process. In the present 
case, thermalization implies that the time-evolved state 
has the same properties as the uniquely defined ther- 
mal equilibrium state whose temperature is such that 
its energy equals the energy of the pump excited system, 
S(t > 4) = Tr[e- H / T « 1! H]/Z. The latter condition defines 
an effective temperature T e g(£ ), and we test for thermal- 
ization by comparing time-dependent expectation val- 
ues of observables 0(t) to thermal expectation values 
0(T cff ) = Tr[e- H / T " !! 0}/Z which arc obtained from an 
independent equilibrium DMFT calculation. To com- 
pute the thermal energy £(T) = Tr[e~ H / T H]/Z which 
is needed for the determination of T e ff, we use the same 
impurity solver as in the nonequilibrium calculation, i.e., 
the same order of the self-consistent hybridization expan- 
sion. 

In Fig. QJ we compare thermal and time-dependent ex- 
pectation values of the double occupancy d(t) = (ni-j-n^), 
which is a purely local quantity and can be measured di- 
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FIG. 2: Relaxation of the double occupancy d(i) after an 
excitation with frequency fi = Un/2 (/3 = 5). Upper panel: 
(a) The difference \d(t) - d(T eS )\, obtained from the OCA 
impurity solver. The pump-amplitudes Eo are given by 3.6701 
(U = 5), 2.4397 (U = 3.5), 2.1325 (U = 3), 1.9123 (U = 2.5), 
1.8198 (U = 2.2), 1.7761 (U = 2), 1.6632 (U = 1.8), 1.2765 
(U — 1.5), such that T e g = 0.5 in all cases. Lower panels: 
Parameters Tth (b) and do (c), obtained from fitting the curves 
in a) for t > 12 with Eq. (|14p (square symbols). Stars and 
triangles show the result of the analogous investigation using 
NCA and the 3rd order hybridization expansion, respectively. 
The time-evolution in 3rd order was restricted to J < 15. 
Black solid lines in b) correspond to the exponential law-22 
Tth = to exp(aU log U), with a = 1.01 and tq = 0.45 (1.2) for 
OCA (NCA). 



rectly in the impurity model. Apparently, a large part of 
the excitation energy goes into the creation of doublon- 
holon pairs, such that d(t) rises by more than 30% dur- 
ing the pump relative to its small value in the initial 
Mott-insulating state. Subsequent to the pump, d{i) 
rapidly settles at a value d s tat 7^ d(T e g) 7 and further re- 
laxation towards the thermal value is hardly visible on 
the scale of Fig. [T] This behavior is similar to that of the 
Hubbard model after an interaction quenchJ^ However, 
a sudden change of the interaction causes pronounced 
2ir/ [/-periodic oscillations which have been used to char- 
acterize the short time behavior in Ref. [3, while these 
oscillations are very weak in the present case, where the 
pump-excitation itself takes already longer than the pe- 
riod 2-k/U. 
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The long-lived nonthermal value of d(t) implies that 
doubly occupied and empty sites do not rapidly re- 
combine in the Hubbard model at large Coulomb 
rcpulsioni 31 i 32 In fact, in experiments with ultra-cold 
atoms it was found that for U 3> W an artificially in- 
creased double occupancy relaxes only on the exponen- 
tially long timcscale r oc exp(a \og(U /W)U /W)^ in 
agreement with theoretical predictions for the lifetime of 
a single doublon in a background of singly occupied sites 
and holes. The reason for this is that multi-particle scat- 
tering events are required to transform a large quantum 
of energy (U) into many kinetic energy excitations of a 
typical amount W, similar to the case of a semiconductor 
where electrons and holes cannot effectively recombine if 
the main relaxation channel is the emission of phonons 
with a much lower energy than the electronic band gap. 
The nonequilibrium DMFT approach now allows us to 
investigate how this decay rate of a single doublon is 
related to the thermalization of a "cloud" of many exci- 
tations, and how the relaxation is modified when U and 
W are comparable in magnitude. 

To answer this question we plot the difference \d(t) — 
d(T c g) | for a series of initial states at various values of the 
Coulomb interaction and temperature T = 0.2 (Fig. [2^). 
The pump frequency is Q = nU/2, and the pump am- 
plitude Eq is tuned such that the effective temperature 
after the pump is given by T e g = 0.5. The time-evolution 
d(t) is initially dominated by a damped transient, which 
subsequently turns into a smooth exponential relaxation 
towards the thermal value d(T e s) . The analysis of this ex- 
ponential tail offers a rigorous way to separate the rapid 
short-time relaxation from the slower long-time thermal- 
ization: From a linear fit to the curves in Fig. [2^ (for 
t > 12), i.e. an ansatz 

d(t) ~ d{T cff ) + d exp(-Vrth), (14) 

we extract both the long-time thermalization rate 1/rth, 
and a value do which measures the "residue" of the ini- 
tial rapid relaxation process. In the limit U ^> W, do 
becomes identical to the apparent nonthermal stationary 
double occupancy d s t&t introduced in the discussion of 
Fig. [T] The numerical effort required to resolve the small 
difference d(t) — d(T e s) restricts the accessible times, such 
that a relaxation over more than one order of magnitude 
can be observed only for a limited range of interactions 
(U = 2, 2.2, 2.5 in Fig. Nevertheless, the behavior of 
Tth and do as a function of U which is obtained on the 
basis of our present results already draws an interesting 
picture for the relaxation of the Hubbard model at strong 
and intermediate coupling, which we will discuss in the 
following two paragraphs. 

For U ^> W, where the total number of excited dou- 
blons is small, we find that the dependence of the ther- 
malization rate r t h on U is indeed consistent with the 
exponential law 3 ^ r t h = To exp(a£/ log U) that was pre- 
dicted for the decay of a single doublon (Fig. EJd). Up 
to numerical accuracy, NCA, OCA, and the third or- 
der of the hybridization expansion give the same large- U 



asymptotics for r t h (a ~ 1) and differ only by the pref- 
actor To- Due to the exponential dependence on [/, the 
thermalization rate is much larger than the inverse band- 
width already for moderately strong interactions. 

When U and W are similar in magnitude, the relax- 
ation changes in two ways: (i) the rate Tth significantly 
decreases to values of the order of a few inverse W and 
starts to deviate from the large-t/ asymptotics, and (ii), 
the residue do extrapolates to zero at a finite value of the 
interaction. While both facts eventually imply thermal- 
ization on a timescale given by the inverse bandwidth, 
they are conceptually quite different. The vanishing of 
the residue do implies the absence of the long-time ex- 
ponential relaxation tail, such that the system follows a 
qualitatively different path to the thermal state. In fact, 
the vanishing of the rapid relaxation residue do provides a 
rigorous way to define a point of "rapid thermalization" , 
as opposed to a regular, two step relaxation. For U < 1.5 
thermalization seems to become slower again, in agree- 
ment with similar findings for the interaction quench to 
the weak-coupling regime , 14 ' 16 but due to the limited nu- 
merical accuracy we cannot fit an exponential tail to the 
curves for U < 1.8 and study this regime systematically. 
That the system approaches a rapid thermalization point 
in the current setup may come as a surprise, because 
neither the sequence of initial states at T = 0.2, nor 
that of the the final states at T = 0.5 crosses a known 
phase-transition of the half-filled Hubbard model. The 
phenomenon is clearly beyond a simple quasi-equilibrium 
explanation, and it is an interesting question whether it 
is related to the "rapid thermalization" at the dynami- 
cal transition in the Hubbard model after an interaction 
quench^ and whether the latter can be characterized in 
a similar way. 



IV. TIME-RESOLVED PHOTOEMISSION 
SPECTROSCOPY 

The results of the previous section raise the question 
how well the different thermalization behaviors at strong 
and intermediate coupling can be distinguished by means 
of spectroscopic methods. In this section we will inves- 
tigate this question for the case of time-resolved pho- 
tocmission spectroscopy. The signal I{u>, t p ) in this case 
is simply defined as the probability for an electron to 
get emitted into an outgoing state with kinetic energy 
hk 2 /2m = ficjprobo — w — $ in response to a probe-pulse 
with frequency cj pro bc at the probe-time t p ($ is the un- 
known work- function). Similar to equilibrium, the pho- 
toemission spectrum can be related to the single-particle 
Greenfunction G < (t 1 if) = ii^^it^Cju^)) on a lattice site 

,■33 

J i 

I(uj,t p ) = -i Jdtdtf S(t)S(t')e iu,( - t '-^G < (t p + t,t p + t'), 

(15) 
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FIG. 3: Time-resolved photoemission spectrum [Gaussian 
probe, 5 = 2, t p = 10] of a pump-excited system (U = 5, 
/3 — 5), using the pump parameters of Fig. [1] The blue dot- 
ted line marked T" e ff = 0.5 is the photoemission spectrum of a 
system at temperature T — 0.5, for the same Gaussian probe 
pulse as for the other two curves. 



where S(t) is the envelope of the probe pulse. For the 
following investigation we use a Gaussian probe S(t) = 
cxp(-£ 2 /2(5 2 )e(3<5 - |t|), where S is the duration. We 
cut off the exponential tails of the probe for \t\ > t c = 
35, such that Eq. (fPo)) can be evaluated for probe times 
t c < t p < i max — t c if the Greenfunction is known for 
0<t,t' <t max . 

For the derivation of Eq. (|T5|) one has to disregard 
matrix element effects and use the sudden approxima- 
tion, which neglects the interaction between electrons 
in the outgoing states and in the solid^ These ap- 
proximations are the same that are commonly made 
in order to compare photoemission spectra in equilib- 
rium to the product of the spectral function A(oS) and 
the Fermi function f(uj). In fact, for an equilibrium 
state, G < (t, t') is just given by the Fourier transform 
G<(t,t') = i J (icjA(w)/(w)e^( t '- t ), such that Eq. flU]) is 
reduced to a convolution of the product A(oj)f(uj) with 
the power density \S(lu)\ 2 = \ J dtS(t)e lult \ 2 of the probe 
envelope S(t), i.e., 7(w) = / du'\S(u)')\ 2 A(u) - u')f(u - 
u> ). Hence A(lu) can in principle be measured with ar- 
bitrary accuracy by choosing a sufficiently long probe- 
pulse. On the other hand, for the noncquilibrium case the 
intrinsic frequency-time uncertainty in Eq. (|15|) makes it 
impossible to measure (or even define) transport of spec- 
tral weight on the timescale of the inverse bandwidth^ 
In the present case, however, we will show that time- 
resolved photoemission can be used to identify the ther- 
malization signatures that have been discussed in the pre- 
vious section. 

In Fig. [3] we reconsider the two pump-experiments of 
Fig. [Don a Mott-insulating initial state (17 = 5, j3 = 5), 
with different pump-frequencies but an equal amount of 
absorbed energy (f2 = ttII/4 and D. = nil/2, T c g = 0.5). 
The photoemission spectrum I(uj,t p ) [Eq. [15])] is broad- 
ened due to a relatively short pulse duration 6 = 2, but 
the Hubbard bands and a gap can still clearly be dis- 
tinguished. Similar to the slow relaxation of the double 
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FIG. 4: Time-resolved photoemission spectrum [Gaussian 
probe, 5 = 2] of a pump-excited system (U = 2.5, j3 = 5), us- 
ing the pump parameters of Fig. [2ji (fi = 57r/4, _Eo=1.9123). 
Main plot: I(u),t p ) for various probe-times t p , compared to 
the spectrum /r e£f (w) of the equilibrium state at T e g = 0.5. 
For better visibility of the upper Hubbard band, data have 
been divided by a Fermi-function with temperature T 1 = 0.55. 
{T' is chosen by hand and is slightly larger than T c a, in 
order to account for the broadening effect due to the finite 
pulse duration.) Inset: Relative difference I(uj,t p ) to JT eff (w) 
as a function of t p , for three probe frequencies wm = 2, 
w (2) = 0.5, W(3) = —0.5. Solid lines show exponential fits 
Jn(w) — lT aff (uj) oc exp(— t/rth), with parameters r-th = 4.70 
(w ( i)), 4.82 (w (2) ), 4.77 (w (3) ). 



occupancy at U = 5, I(uj, t p ) becomes almost indepen- 
dent of t p after the decay of the initial transient for all 
times accessible with the OCA impurity solver. In the 
figure, we compare the spectra of the pump-excited sys- 
tem at tp = 10 to equilibrium spectra both at the initial 
temperature (T = 0.2) and at T = T c g. The relative 
differences between these curves are most pronounced in 
the upper Hubbard band (UHB). Due to thermal fluctu- 
ations, the latter acquires much more weight for T = T e g 
than for the initial equilibrium state T = 0.2. For the 
two pump-excited states, the weight of the UHB is com- 
parable to the weight of the UHB at T = T e g, but its 
distribution depends on the pumping process in such a 
way that the spectral density at high energies is increased 
for the high-frequency pump. 

The differences between the curves in Fig. [3] clearly 
confirm the nonthermal nature of the excited states, 
which had already been deduced from the value of the 
double occupancy. However, a thermal state at T = T c g 
might not be accessible in a typical pump-probe exper- 
iment, because electronic temperatures can be so high 
that heating the whole sample to those temperatures 
would result in its destruction, or at least in a strong 
modification of its properties. It is thus important to note 
that no quantitative comparison to the thermal state at 
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T = T c ff is required in order to identify a lack of thermal- 
ization of the electrons, but it is enough to demonstrate 
that a state still has detailed memory on the way it was 
prepared. In the present case this memory becomes ev- 
ident from the fact that the spectral weight distribution 
depends on the pump-frequency, for systems which have 
absorbed the same amount of energy. 

At intermediate coupling, the time-evolution of the 
double occupancy has indicated a more rapid thermaliza- 
tion of the system, and we will now demonstrate that this 
result is confirmed by the behavior of the spectral func- 
tion. For this purpose we choose a state in the crossover 
regime (U = 2.5, /3 = 5), excite it by a pulse of frequency 
ft = Un/2 to the energy corresponding to T e g = 0.5, and 
compare the spectrum I(iv,t p ) to the spectrum /r eff (a;) 
of the thermal equilibrium state at T = T e g (Fig. @|. At 
these parameters, the Mott gap is already closed, but the 
Hubbard bands are still separated by a dip in the spectral 
function. In the photoemission spectrum, the UHB ap- 
pears merely as a shoulder on the upper edge, but it can 
be enhanced by dividing the spectrum by a Fermi func- 
tion. For small probe-times, I(ui,t p ) differs from It b{{ (oj) 
by a shift of the UHB to larger frequencies, but in con- 
trast to the results for U = 5 we now observe a relaxation 
of the spectrum towards the thermal equilibrium spec- 
trum at later t p . In the inset of Fig. we plot the relative 
difference AJ = t p ) — ^r e£f (w)]//T e ff ( w ) as a function 
of probe time t p for three frequencies, cJm = 2 (above 
the center of the UHB), cjt 2 ) — 0-5 (below the center of 
the UHB), and W( 2 ) = —0.5 . Due to the shift of spec- 
tral weight from higher to lower frequencies, I(to^,t p ) 
decreases with t p , while /(w(2),i P ) increases. A fit of the 
curves with an exponential decay gives relaxation times 
Ttj, = 4.70 for W(i), r t h = 4.82 for W( 2 ), and r t h = 4.77 for 
W(3), comparable to the relaxation time extracted from 
d(t) for the same parameters, r t h=4.87 (Fig. [2]). 

V. CONCLUSION 

In conclusion, we have investigated the relaxation of 
a Mott insulator subsequent to a strong laser excitation 
over the Hubbard gap. After an initial strongly damped 
transient, the time-evolution of the double occupancy 
turns into an exponential relaxation towards its thermal- 
ized value at large times. For U 3> W, the thermalization 
time grows rapidly with U/W, consistent with the expo- 
nentially small decay rate for one doublon in front of a 
background of singly occupied sites and holes^ When U 
is comparable to W, on the other hand, the double occu- 
pancy thermalizes within a few times the inverse band- 
width. This more efficient thermalization does not only 
become manifest through a decrease of the thermaliza- 
tion rate 1/ith, but moreover, the system approaches a 
point at which the exponential relaxation is entirely ab- 
sent, and the thermal state is approached along a quali- 
tatively different, "rapid" pathway. 



Exponential thermalization arises when the relevant 
kinetic equations can be linearized around the thermal 
state, i.e., when relaxation close to a thermal equilibrium 
state is described in terms of small deviations from that 
state. The fact that the long-time relaxation of the dou- 
ble occupancy can be understood from the decay rate of 
single doublons fits into this scheme. On the other hand, 
the appearance of the rapid thermalization instead of the 
exponential relaxation is quite remarkable, in particular 
as it appears at a value of U that is far from any known 
phase transition in equilibrium. It remains to be seen 
whether similar results can be obtained for weak cou- 
pling (which is not accessible with the impurity solver 
used in this work), where the long-time behavior should 
be described by a Boltzmann equation^ and whether it 
is related to the "rapid thermalization" at the dynami- 
cal transition after an interaction quench in the Hubbard 
model. In any case, the analysis of the long-time tail in 
terms of the "residue" do used in this paper allows to rig- 
orously define a "rapid thermalization" , and thus should 
be useful for a future analysis of the dynamical transition 
after the interaction quench in the Hubbard model and 
other models. 

In this work we have also demonstrated that the dif- 
ferences in thermalization behavior are in principle ob- 
servable by means of time-resolved photoemission spec- 
troscopy. Quantitative differences between the time- 
dependent spectrum and equilibrium spectra are small, 
but nonthermal behavior of the pump-excited system 
may also be detected from qualitative features, such as 
the memory of the electronic system on the precise form 
of the excitation. We suppose that similar information 
can be drawn from optical measurements, since both the 
conductivity and the photoemission spectrum are derived 
from the same single-particle Green function (up to ver- 
tex corrections which can be controlled to some extent by 
varying the relative polarization of pump and probe) ^ 
Of course, most materials are more complicated than the 
simple one-band Hubbard model studied in this work. 
But our results suggest that a similar pronounced de- 
pendence of the relaxation on the interaction parameter 
occurs in more complex models. For example, it would be 
interesting to see whether the relaxation of photo-doped 
metallic phases to insulating phases in charge transfer 
insulators^ is subject to a similar strong dependence on 
the system parameters as the thermalization of doublons 
in the half-filled single-band system. 
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